## Project: Autocracy UPR project ----------------------------------------------
##        05 Appendix             ----------------------------------------------
## Author:
##   - Chun-Young Park (UGA)
##   - Sang-Hoon Park (UofSC)
##
## Updated: Aug 31 2024
if (!require("devtools")) install.packages("devtools")
if (!require("pacman")) install.packages("pacman")
if (!require("rqog")) remotes::install_github("ropengov/rqog")
pacman::p_load(rqog, rio, tidyverse, dendextend, splitstackshape, psych, tidyverse)
if (!require("dendextendRcpp")) devtools::install_github('talgalili/dendextendRcpp')
if (!require("vdemdata")) devtools::install_github("vdeminstitute/vdemdata")
if (!require("ggrepel")) install.packages("ggrepel")
if (!require("patchwork")) install.packages("patchwork")
if (!require("futurevisions")) devtools::install_github("JoeyStanley/futurevisions")
if (!require("ERT")) devtools::install_github("vdeminstitute/ERT")
if (!require("wesanderson")) devtools::install_github("karthik/wesanderson")

showtext::showtext_auto(T)

## theme_set
theme_clean <- function() {
  theme_minimal(base_family = "Barlow Semi Condensed") +
    theme(panel.grid.minor = element_blank(),
          plot.title = element_text(face = "bold", color = "black"),
          plot.subtitle = element_text(face = "bold", color = "black"),
          axis.title = element_text(family = "Barlow Semi Condensed Medium", color = "black"),
          axis.text = element_text(color = "black"),
          strip.text = element_text(family = "Barlow Semi Condensed",
                                    face = "bold", size = rel(1), hjust = 0,
                                    color = "black"),
          strip.background = element_rect(fill = "white", color = NA),
          plot.caption = element_text(hjust = 0, color = "black"),
          legend.position = "bottom")
}

ggplot2::theme_set(
  theme_clean()
)

sample <- readRDS("data/data for analysis/sample.RDS")
# If you do not load RDS file, you have to run codes in 05_analysis.R

ologit_model1 <- readRDS("data/estimation/ologit_model1.RDS");
ologit_model2 <- readRDS("data/estimation/ologit_model2.RDS");
ologit_model3 <- readRDS("data/estimation/ologit_model3.RDS");
ologit_model4 <- readRDS("data/estimation/ologit_model4.RDS")

## A. Descriptive Statistics ---------------------------------------------------
### Table A.1: Descriptive statistics of sample --------------------------------
sample |> dplyr::select(
  Specificity = severity, 
  `Democratic Dyads` = dd_lexi, 
  `Autocratic Dyads` = aa_lexi, 
  `Auto$_\\text{Reviewer}$-Demo$_\\text{SUR}$ Mixed Dyads` = ad_lexi, 
  `Demo$_\\text{Reviewer}$-Auto$_\\text{SUR}$ Mixed Dyads` = da_lexi,
  Year = year,
  `EDI difference$_\\text{Reviewer - SUR}$` = politysenderminustarget, 
  `Geopolitical affinity` = idealdistance, 
  Alliance = alliance, 
  `Ln(Combined GDPpc$_\\mathrm{SuR-Reviewer}$)` = share_combined, 
  `Reviewer reviewed`= reviewer_revieweed, 
  `HRC Member$_\\mathrm{SuR}$` = to_mem, 
  `HRC Member$_\\mathrm{Reviewer}$` = from_mem, 
  `Same region` = regionalism, 
  `Human rights score$_\\mathrm{SUR}$` = target_hrs, 
  `Human rights score$_\\mathrm{Reviewer}$` = sender_hrs, 
  Migration = migrantscat, 
  `Socio-economic rights` = socioecon, 
  `Vulnerable populations` = vulnerable, 
  `Women, children tracking` = womenchild, 
  `Physical integrity rights` = physint, 
  `Justice` = justicecat, 
  `Speech political participation` = political
) |> drop_na() |> psych::describe() |> 
  as.data.frame() |> rownames_to_column() |> 
  dplyr::select(
    Variables = rowname,
    `No. of obs` = n,
    Mean = mean,
    Median = median,
    `Std.dev.` = sd,
    Min = min, Max = max
  ) |> kableExtra::kbl("latex", 
                    booktabs = T, format.args = list(big.mark = ","),
                    longtable = F,
                    digits = 1) |> 
  kableExtra::kable_paper(font_size = 9) |> 
  kableExtra::save_kable(format = "latex",
                         "tables/appendix_table_a1.tex")

### Table A.2: Descriptive statistics of sample with 1-year lagged bilateral trades -------
sample |> dplyr::select(
  Specificity = severity, 
  `Democratic Dyads` = dd_lexi, 
  `Autocratic Dyads` = aa_lexi, 
  `Auto$_\\text{Reviewer}$-Demo$_\\text{SUR}$ Mixed Dyads` = ad_lexi, 
  `Demo$_\\text{Reviewer}$-Auto$_\\text{SUR}$ Mixed Dyads` = da_lexi,
  Year = year,
  `EDI difference$_\\text{Reviewer - SUR}$` = politysenderminustarget, 
  `Geopolitical affinity` = idealdistance, 
  Alliance = alliance, 
  `L.ln(trade$_{t-1}$)` = lagtrade1, 
  `Reviewer reviewed`= reviewer_revieweed, 
  `HRC Member$_\\mathrm{SuR}$` = to_mem, 
  `HRC Member$_\\mathrm{Reviewer}$` = from_mem, 
  `Same region` = regionalism, 
  `Human rights score$_\\mathrm{SUR}$` = target_hrs, 
  `Human rights score$_\\mathrm{Reviewer}$` = sender_hrs, 
  Migration = migrantscat, 
  `Socio-economic rights` = socioecon, 
  `Vulnerable populations` = vulnerable, 
  `Women, children tracking` = womenchild, 
  `Physical integrity rights` = physint, 
  `Justice` = justicecat, 
  `Speech political participation` = political
) |> drop_na() |> psych::describe() |> 
  as.data.frame() |> rownames_to_column() |> 
  dplyr::select(
    Variables = rowname,
    `No. of obs` = n,
    Mean = mean,
    Median = median,
    `Std.dev.` = sd,
    Min = min, Max = max
  ) |> kableExtra::kbl("latex", 
                       booktabs = T, format.args = list(big.mark = ","),
                       longtable = F,
                       digits = 1) |> 
  kableExtra::kable_paper(font_size = 9) |> 
  kableExtra::save_kable(format = "latex",
                         "tables/appendix_table_a2.tex")

### Table A.3: Descriptive statistics of sample with 3-year lagged bilateral trades -------
sample |> dplyr::select(
  Specificity = severity, 
  `Democratic Dyads` = dd_lexi, 
  `Autocratic Dyads` = aa_lexi, 
  `Auto$_\\text{Reviewer}$-Demo$_\\text{SUR}$ Mixed Dyads` = ad_lexi, 
  `Demo$_\\text{Reviewer}$-Auto$_\\text{SUR}$ Mixed Dyads` = da_lexi,
  Year = year,
  `EDI difference$_\\text{Reviewer - SUR}$` = politysenderminustarget, 
  `Geopolitical affinity` = idealdistance, 
  Alliance = alliance, 
  `L.ln(trade$_{t-3}$)` = lagtrade3, 
  `Reviewer reviewed`= reviewer_revieweed, 
  `HRC Member$_\\mathrm{SuR}$` = to_mem, 
  `HRC Member$_\\mathrm{Reviewer}$` = from_mem, 
  `Same region` = regionalism, 
  `Human rights score$_\\mathrm{SUR}$` = target_hrs, 
  `Human rights score$_\\mathrm{Reviewer}$` = sender_hrs, 
  Migration = migrantscat, 
  `Socio-economic rights` = socioecon, 
  `Vulnerable populations` = vulnerable, 
  `Women, children tracking` = womenchild, 
  `Physical integrity rights` = physint, 
  `Justice` = justicecat, 
  `Speech political participation` = political
) |> drop_na() |> psych::describe() |> 
  as.data.frame() |> rownames_to_column() |> 
  dplyr::select(
    Variables = rowname,
    `No. of obs` = n,
    Mean = mean,
    Median = median,
    `Std.dev.` = sd,
    Min = min, Max = max
  ) |> kableExtra::kbl("latex", 
                       booktabs = T, format.args = list(big.mark = ","),
                       longtable = F,
                       digits = 1) |> 
  kableExtra::kable_paper(font_size = 9) |> 
  kableExtra::save_kable(format = "latex",
                         "tables/appendix_table_a3.tex")

## B. Model Specifications (Benchmark and Full Models) ---------------------------------------------
### Table B.1 shows the models of specifications with no controls, the share of combined GDP -------

texreg::texreg(list( 
  ologit_model1, ologit_model2,
  ologit_model3, ologit_model4),
  omit.coef = "as.factor",
  custom.model.names=c(
    "Model 1",
    "Model 2", 
    "Model 3",
    "Model 4"),
  custom.coef.map = list(
    "aa_lexi" = "Autocracy Dyads",
    "da_lexi" = "Demo$_{\\mathrm{Reviewer}}$-Auto$_{\\mathrm{SuR}}$ Mixed Dyads", 
    "ad_lexi" = "Auto$_{\\mathrm{Reviewer}}$-Demo(SUR)",
    "1|2" = "(Cut 1)",
    "2|3" = "(Cut 2)",
    "politysenderminustarget" = "Joint Democracy", 
    "share_combined" = "Share of Combined GDP",
    "alliance" = "Alliance", 
    "idealdistance" =  "Geopolitical Affinity",
    "reviewer_revieweed" = "Reviewed Reviewer",
    "to_mem" = "HRC Member$_{\\mathrm{SuR}}$", 
    "from_mem" = "HRC Member$_{\\mathrm{Reviewer}}$",
    "regionalism"= "Same Region",
    "sender_hrs" = "Human Rights Score$_{\\mathrm{Reviewer}}$",
    "target_hrs" = "Human Rights Score$_{\\mathrm{SuR}}$",
    "migrantscat"= "Migration", 
    "socioecon" = "Socio-Economic Rights", 
    "vulnerable" = "Vulnerable Populations",
    "womenchild" =   "Women, Children \\& Trafficking", 
    "physint" = "Physical Integrity Rights",
    "justicecat" =   "Justice", 
    "political" = "Speech \\& Political Participation",
    "lagtrade1" =   "L.ln(trade$_{t−1}$)",
    "lagtrade3" =   "L.ln(trade$_{t−3}$)"),
  reorder.coef = c(1, 2, 3, 7, 23, 24, 6, 8, 9, 11, 10, 12, 13, 14,
                   15, 16, 17, 18, 19, 20, 21, 22, 4, 5),
  single.row = TRUE,
  include.rsquared=FALSE, 
  custom.gof.rows=list(`Country-fixed` = c("Yes", "Yes", "Yes", "Yes"),
                       `Year-fixed` = c("Yes", "Yes", "Yes", "Yes")),
  stars=c(0.001, 0.01, 0.05), 
  #custom.note="%stars. Year effects not significant.",
  caption=
    "\\label{tab1} Ordered logit models on UPR recommendations shaming behaviors",
  file = "tables/appendix_table_b1.tex")

## C. Alternative Measures of Political Regimes ---------------------------------------------
vdemdata::vdem -> vdem
vdem |> 
  mutate(
    polyarchy = if_else(
      is.na(v2x_polyarchy), NA_integer_, if_else(
        v2x_polyarchy < 0.5, 0L, 1L)),
    vdem = if_else(
      is.na(v2x_regime), NA_integer_, if_else(
        v2x_regime %in% c(0, 1), 0L, 1L)),
    brm = e_boix_regime,
    lied5 = if_else(
      is.na(e_lexical_index), NA_integer_, if_else(
        e_lexical_index < 6, 0L, 1L)),
    lied4 = if_else(
      is.na(e_lexical_index), NA_integer_, if_else(
        e_lexical_index < 5, 0L, 1L)),
    lied3 = if_else(
      is.na(e_lexical_index), NA_integer_, if_else(
        e_lexical_index < 4, 0L, 1L)),
    polity = if_else(
      is.na(e_polity2), NA_integer_, if_else(
        e_polity2 < 7 , 0L, 1L))
  ) |> dplyr::select(year, polyarchy,vdem, brm, lied5, lied4, lied3, polity) |> 
  pivot_longer(
    cols = c("polyarchy":"polity"),
    names_to = "measurements",
    values_to = "scores"
  ) |> dplyr::filter(year > 2007) -> 
  append_measurements_regime

append_measurements_regime |> 
  group_by(year, measurements, scores) |> count() |> drop_na() |> 
  mutate(regime = factor(scores, levels = c(0, 1),
                         labels = c("Autocracy", "Democracy"))) |> 
  ungroup() |> 
  mutate(n = if_else(n < 2, NA_integer_, n)) -> append_measurements_regime

data_ends <- append_measurements_regime |> 
  mutate(measurements = 
           factor(measurements,
                  levels = c("brm", "polyarchy", "vdem", "lied3", "lied4", "lied5", "polity"),
                  labels = c("Boix et al. Dichotomous",
                             "V-Dem EDI > 0.5", 
                             "RoW: Closed/electoral",
                             "LIED > 3", "LIED > 4", "LIED > 5",
                             "POLITY2 > 6"))) |> 
  group_by(measurements, regime) |> drop_na(scores, n, regime) |> 
  dplyr::filter(year %in% max(year)) |> ungroup() 


append_measurements_regime |> ungroup() |>
  mutate(measurements = 
           factor(measurements,
                  levels = c("brm", "polyarchy", "vdem",  "lied3", "lied4", "lied5", "polity"),
                  labels = c("Boix et al. Dichotomous",
                             "V-Dem EDI > 0.5", 
                             "RoW: Closed/electoral",
                             "LIED > 3", "LIED > 4", "LIED > 5",
                             "POLITY2 > 6"))) |> 
  ggplot(aes(x = year, y = n, color = measurements)) + 
  geom_line(show.legend = F,
            position = position_dodge2(1)) +
  geom_point(fill = "white", size = 3, show.legend = F,
             position = position_dodge2(1)) +
  facet_wrap(~regime) + 
  geom_label_repel(
    aes(label = measurements), data = data_ends,
    fontface ="plain", size = 3, show.legend = F) +
  scale_x_continuous(breaks = c(seq(2008, 2023, 1)),
                     labels = c("08", "09", "10", "11", "12", "13", "14", "15",
                                "16", "17", "18", "19", "20", "21", "22", "23")) +
  ggsci::scale_color_lancet() +
  labs(y = "Frequencies\n", x = "\nYear")

#### Save the patterns of regime classification for different indicators
### Figure C.1: Annual change in the numbers of different political regime by ----
### different measurements
ggsave("figures/appendix_figure_c1.pdf", width = 12, height = 6, dpi = "retina")

## D. Bivariate Relationship between Regime Dyads and UPR severity ------------
### Figure D.1: The proportion of recommendation severity by different dyads ----
###             based on different regime meaSURements
#### LIED < 5; Democracy and Dictatorship
sample |> drop_na(severity) |>
  drop_na(regime_dyads_lexi4) |>
  group_by(regime_dyads_lexi4, severity) |> 
  dplyr::select(7) |> summarize(total = n()) |> ungroup() |>
  group_by(regime_dyads_lexi4) |> 
  mutate(regime_dyads_lexi4 = 
           factor(regime_dyads_lexi4,
                  levels = c("Democracy Dyads", "Autocracy Dyads",
                             "Auto(REV)-Demo(SUR) Mixed Dyads", 
                             "Demo(REV)-Auto(SUR) Mixed Dyads"))) |> 
  mutate(sum = sum(total, na.rm = T),
         ratio = total/sum,
         severity_fa = case_when(
           severity == 1L ~ "Praising",
           severity == 2L ~ "Neutral",
           severity == 3L ~ "Shaming",
           T ~ NA_character_),
         severity_fa = factor(severity_fa, levels = c(
           "Praising", "Neutral", "Shaming"))) |>
  ggplot(aes(x = regime_dyads_lexi4, fill = severity_fa, y = ratio)) + 
  geom_bar(stat = "identity", position=position_dodge(0.9)) + 
  labs(x = "", y = NULL,
       subtitle = "Autocracy: LIED of Skaaning et al. (2015) < 4") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0, 0.6)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_fill_manual(values = c(futurevisions("mars")[3],
                               futurevisions("mars")[2],
                               futurevisions("mars")[1])) + 
  theme(legend.position = "bottom", legend.title = element_blank()) ->
  p1

##### LIED < 5; Democracy and Dictatorship
sample |> drop_na(severity) |>
  drop_na(regime_dyads_lexi5) |>
  group_by(regime_dyads_lexi5, severity) |> 
  dplyr::select(7) |> summarize(total = n()) |> ungroup() |>
  group_by(regime_dyads_lexi5) |> 
  mutate(regime_dyads_lexi5 = 
           factor(regime_dyads_lexi5,
                  levels = c("Democracy Dyads", "Autocracy Dyads",
                             "Auto(REV)-Demo(SUR) Mixed Dyads", 
                             "Demo(REV)-Auto(SUR) Mixed Dyads"))) |> 
  mutate(sum = sum(total, na.rm = T),
         ratio = total/sum,
         severity_fa = case_when(
           severity == 1L ~ "Praising",
           severity == 2L ~ "Neutral",
           severity == 3L ~ "Shaming",
           T ~ NA_character_),
         severity_fa = factor(severity_fa, levels = c(
           "Praising", "Neutral", "Shaming"))) |>
  ggplot(aes(x = regime_dyads_lexi5, fill = severity_fa, y = ratio)) + 
  geom_bar(stat = "identity", position=position_dodge(0.9)) + 
  labs(x = "", y = NULL,
       subtitle = "Autocracy: LIED of Skaaning et al. (2015) < 5") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0, 0.6)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_fill_manual(values = c(futurevisions("mars")[3],
                               futurevisions("mars")[2],
                               futurevisions("mars")[1])) + 
  theme(legend.position = "bottom", legend.title = element_blank()) ->
  p2

##### Polyarchy < 0.5; Democracy and Dictatorship
sample |> drop_na(severity) |>
  drop_na(regime_dyads) |>
  group_by(regime_dyads, severity) |> 
  dplyr::select(7) |> summarize(total = n()) |> ungroup() |>
  group_by(regime_dyads) |> 
  mutate(regime_dyads = 
           factor(regime_dyads,
                  levels = c("Democracy Dyads", "Autocracy Dyads",
                             "Auto(REV)-Demo(SUR) Mixed Dyads", 
                             "Demo(REV)-Auto(SUR) Mixed Dyads"))) |> 
  mutate(sum = sum(total, na.rm = T),
         ratio = total/sum,
         severity_fa = case_when(
           severity == 1L ~ "Praising",
           severity == 2L ~ "Neutral",
           severity == 3L ~ "Shaming",
           T ~ NA_character_),
         severity_fa = factor(severity_fa, levels = c(
           "Praising", "Neutral", "Shaming"))) |>
  ggplot(aes(x = regime_dyads, fill = severity_fa, y = ratio)) + 
  geom_bar(stat = "identity", position=position_dodge(0.9)) + 
  labs(x = "", y = NULL,
       subtitle = "Autocracy: V-Dem's Electoral Democracy Index < 0.5") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0, 0.6)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_fill_manual(values = c(futurevisions("mars")[3],
                               futurevisions("mars")[2],
                               futurevisions("mars")[1])) + 
  theme(legend.position = "bottom", legend.title = element_blank()) ->
  p3

##### POLITY4 < 7; Democracy and Dictatorship
sample |> drop_na(severity) |>
  drop_na(regime_dyads_pol4) |>
  group_by(regime_dyads_pol4, severity) |> 
  dplyr::select(7) |> summarize(total = n()) |> 
  ungroup() |>
  group_by(regime_dyads_pol4) |> 
  mutate(regime_dyads_pol4 = 
           factor(regime_dyads_pol4,
                  levels = c("Democracy Dyads", "Autocracy Dyads",
                             "Auto(REV)-Demo(SUR) Mixed Dyads", 
                             "Demo(REV)-Auto(SUR) Mixed Dyads"))) |> 
  mutate(sum = sum(total, na.rm = T),
         ratio = total/sum,
         severity_fa = case_when(
           severity == 1L ~ "Praising",
           severity == 2L ~ "Neutral",
           severity == 3L ~ "Shaming",
           T ~ NA_character_),
         severity_fa = factor(severity_fa, levels = c(
           "Praising", "Neutral", "Shaming"))) |>
  ggplot(aes(x = regime_dyads_pol4, 
             fill = severity_fa, y = ratio)) + 
  geom_bar(stat = "identity", position=position_dodge(0.9)) + 
  labs(x = "", y = NULL,
       subtitle = "Autocracy: POLITY2 < 7") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0, 0.6)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_fill_manual(values = c(futurevisions("mars")[3],
                               futurevisions("mars")[2],
                               futurevisions("mars")[1])) + 
  theme(legend.position = "bottom", legend.title = element_blank()) ->
  p4

#### Boix, Rosato, and Miller; Democracy and Dictatorship
sample |> drop_na(severity) |>
  drop_na(regime_dyads_brm) |>
  group_by(regime_dyads_brm, severity) |> 
  dplyr::select(7) |> summarize(total = n()) |> ungroup() |>
  group_by(regime_dyads_brm) |> 
  mutate(regime_dyads_brm = 
           factor(regime_dyads_brm,
                  levels = c("Democracy Dyads", "Autocracy Dyads",
                             "Auto(REV)-Demo(SUR) Mixed Dyads", 
                             "Demo(REV)-Auto(SUR) Mixed Dyads"))) |> 
  mutate(sum = sum(total, na.rm = T),
         ratio = total/sum,
         severity_fa = case_when(
           severity == 1L ~ "Praising",
           severity == 2L ~ "Neutral",
           severity == 3L ~ "Shaming",
           T ~ NA_character_),
         severity_fa = factor(severity_fa, levels = c(
           "Praising", "Neutral", "Shaming"))) |>
  ggplot(aes(x = regime_dyads_brm, fill = severity_fa, y = ratio)) + 
  geom_bar(stat = "identity", position=position_dodge(0.9)) + 
  labs(x = "", y = NULL,
       subtitle = "Autocracy: Dichotomous meaSURe by Boix et al. (2012).") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0, 0.6)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
  scale_fill_manual(values = c(futurevisions("mars")[3],
                               futurevisions("mars")[2],
                               futurevisions("mars")[1])) + 
  theme(legend.position = "bottom", legend.title = element_blank()) ->
  p5

#### V-Dem; Democracy and Dictatorship
sample |> drop_na(severity) |>
  drop_na(regime_dyads_vdem) |>
  group_by(regime_dyads_vdem, severity) |> 
  dplyr::select(7) |> summarize(total = n()) |> 
  ungroup() |>
  group_by(regime_dyads_vdem) |> 
  mutate(regime_dyads_vdem = 
           factor(regime_dyads_vdem,
                  levels = c("Democracy Dyads", "Autocracy Dyads",
                             "Auto(REV)-Demo(SUR) Mixed Dyads", 
                             "Demo(REV)-Auto(SUR) Mixed Dyads"))) |> 
  mutate(sum = sum(total, na.rm = T),
         ratio = total/sum,
         severity_fa = case_when(
           severity == 1L ~ "Praising",
           severity == 2L ~ "Neutral",
           severity == 3L ~ "Shaming",
           T ~ NA_character_),
         severity_fa = factor(severity_fa, levels = c(
           "Praising", "Neutral", "Shaming"))) |>
  ggplot(aes(x = regime_dyads_vdem, 
             fill = severity_fa, y = ratio)) + 
  geom_bar(stat = "identity", position=position_dodge(0.9)) + 
  labs(x = "", y = NULL,
       subtitle = "Autocracy: Closed and electoral Autocracy") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0, 0.6)) + 
  scale_x_discrete(labels = function(x) str_wrap(x, width = 10)) +
scale_fill_manual(values = c(futurevisions::futurevisions("mars")[3], 
                             futurevisions::futurevisions("mars")[2], 
                             futurevisions::futurevisions("mars")[1])) +
  theme(legend.position = "bottom", legend.title = element_blank()) ->
  p6


p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2, guides = "collect") +
  patchwork::plot_annotation(
    caption = str_wrap("Note: The y-axis represents the percentage of comments within each dyad type, and the x-axis categorizes the different dyad types. The color-coded bars show the proportion of comments classified as Praising (blue), Neutral (yellow), and Shaming (red). 
                       In the mixed dyads, the first label indicates the regime type of the reviewer, while the second label denotes the regime type of the state under review. For example, 'Auto-Demo Mixed Dyads' means that an autocracy is reviewing a democracy.", 
                       125, exdent = 10)) &
  theme(
    legend.title = element_blank(),
    legend.text = element_text(size = rel(1.3)),
    legend.margin=margin(0,0,0,0),
    legend.box.margin=margin(0, 0, 20, 0),
    plot.caption = element_text(size = rel(1.3)))

ggsave("figures/appendix_figure_d1.pdf", width = 12, height = 10, dpi = "retina")

## E. Robustness Check ---------------------------------------------------------
### E.1. Alternative Measurements for Regime Types -----------------------------
sample |> dplyr::select(
  severity,   
  aa_poly, da_poly, ad_poly, 
  politysenderminustarget, 
  lagtrade1, lagtrade3,
  ln_combined_gdppc, share_combined,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3, lagmeantrade,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_poly
gc()

ologit_model_poly0 <- MASS::polr(
  as.factor(severity) ~ 
    aa_poly + da_poly + ad_poly +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_poly ,Hess = T,
  method = "logistic")

ologit_model_poly0 <- MASS::polr(
  as.factor(severity) ~ 
    aa_poly + da_poly + ad_poly +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_poly, Hess = T,
  method = "logistic")
gc()

ologit_model_poly1 <- MASS::polr(
  as.factor(severity) ~ 
    aa_poly + da_poly + ad_poly +
    politysenderminustarget + 
    share_combined +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_poly, Hess = T,
  method = "logistic")
gc()

ologit_model_poly2 <- MASS::polr(
  as.factor(severity) ~ 
    aa_poly + da_poly + ad_poly +
    politysenderminustarget + 
    lagtrade1 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_poly, Hess = T,
  method = "logistic")
gc()

ologit_model_poly3 <- MASS::polr(
  as.factor(severity) ~ 
    aa_poly + da_poly + ad_poly +
    politysenderminustarget + 
    lagtrade3 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_poly, Hess = T,
  method = "logistic")
gc()

##### RoW 
sample |> dplyr::select(
  severity,   
  aa_vdem, da_vdem, ad_vdem, 
  politysenderminustarget, 
  lagtrade1, lagtrade3,
  ln_combined_gdppc, share_combined,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3, lagmeantrade,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_vdem
gc()

ologit_model_vdem0 <- MASS::polr(
  as.factor(severity) ~ 
    aa_vdem + da_vdem + ad_vdem +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_vdem ,Hess = T,
  method = "logistic")

ologit_model_vdem1 <- MASS::polr(
  as.factor(severity) ~ 
    aa_vdem + da_vdem + ad_vdem +
    politysenderminustarget + 
    share_combined +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_vdem, Hess = T,
  method = "logistic")
gc()

ologit_model_vdem2 <- MASS::polr(
  as.factor(severity) ~ 
    aa_vdem + da_vdem + ad_vdem +
    politysenderminustarget + 
    lagtrade1 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_vdem, Hess = T,
  method = "logistic")
gc()

ologit_model_vdem3 <- MASS::polr(
  as.factor(severity) ~ 
    aa_vdem + da_vdem + ad_vdem +
    politysenderminustarget + 
    lagtrade3 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_vdem, Hess = T,
  method = "logistic")
gc()

##### POLITY IV 
sample |> dplyr::select(
  severity,   
  aa_pol4, da_pol4, ad_pol4, 
  politysenderminustarget, 
  lagtrade1, lagtrade3,
  ln_combined_gdppc, share_combined,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3, lagmeantrade,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_pol4
gc()

ologit_model_pol40 <- MASS::polr(
  as.factor(severity) ~ 
    aa_pol4 + da_pol4 + ad_pol4 +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_pol4 ,Hess = T,
  method = "logistic")
gc()

ologit_model_pol41 <- MASS::polr(
  as.factor(severity) ~ 
    aa_pol4 + da_pol4 + ad_pol4 +
    politysenderminustarget + 
    share_combined +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_pol4, Hess = T,
  method = "logistic")
gc()

ologit_model_pol42 <- MASS::polr(
  as.factor(severity) ~ 
    aa_pol4 + da_pol4 + ad_pol4 +
    politysenderminustarget + 
    lagtrade1 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_pol4, Hess = T,
  method = "logistic")
gc()

ologit_model_pol43 <- MASS::polr(
  as.factor(severity) ~ 
    aa_pol4 + da_pol4 + ad_pol4 +
    politysenderminustarget + 
    lagtrade3 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_pol4, Hess = T,
  method = "logistic")
gc()

##### BRM 
sample |> dplyr::select(
  severity,   
  aa_brm, da_brm, ad_brm, 
  politysenderminustarget, 
  lagtrade1, lagtrade3,
  ln_combined_gdppc, share_combined,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3, lagmeantrade,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_brm
gc()

ologit_model_brm0 <- MASS::polr(
  as.factor(severity) ~ 
    aa_brm + da_brm + ad_brm +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_brm ,Hess = T,
  method = "logistic")
gc()

ologit_model_brm1 <- MASS::polr(
  as.factor(severity) ~ 
    aa_brm + da_brm + ad_brm +
    politysenderminustarget + 
    share_combined +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_brm, Hess = T,
  method = "logistic")
gc()

ologit_model_brm2 <- MASS::polr(
  as.factor(severity) ~ 
    aa_brm + da_brm + ad_brm +
    politysenderminustarget + 
    lagtrade1 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_brm, Hess = T,
  method = "logistic")
gc()

ologit_model_brm3 <- MASS::polr(
  as.factor(severity) ~ 
    aa_brm + da_brm + ad_brm +
    politysenderminustarget + 
    lagtrade3 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_brm, Hess = T,
  method = "logistic")
gc()

##### LIED < 5 
sample |> dplyr::select(
  severity,   
  aa_lied5, da_lied5, ad_lied5, 
  politysenderminustarget, 
  lagtrade1, lagtrade3,
  ln_combined_gdppc, share_combined,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3, lagmeantrade,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_lied5
gc()

ologit_model_lied50 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied5 + da_lied5 + ad_lied5 +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied5 ,Hess = T,
  method = "logistic")
gc()

ologit_model_lied51 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied5 + da_lied5 + ad_lied5 +
    politysenderminustarget + 
    share_combined +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied5, Hess = T,
  method = "logistic")
gc()

ologit_model_lied52 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied5 + da_lied5 + ad_lied5 +
    politysenderminustarget + 
    lagtrade1 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied5, Hess = T,
  method = "logistic")
gc()

ologit_model_lied53 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied5 + da_lied5 + ad_lied5 +
    politysenderminustarget + 
    lagtrade3 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied5, Hess = T,
  method = "logistic")
gc()

##### LIED < 4 
sample |> dplyr::select(
  severity,   
  aa_lied4, da_lied4, ad_lied4, 
  politysenderminustarget, 
  lagtrade1, lagtrade3,
  ln_combined_gdppc, share_combined,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3, lagmeantrade,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_lied4
gc()

ologit_model_lied40 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied4 + da_lied4 + ad_lied4 +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied4 ,Hess = T,
  method = "logistic")
gc()

ologit_model_lied41 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied4 + da_lied4 + ad_lied4 +
    politysenderminustarget + 
    share_combined +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied4, Hess = T,
  method = "logistic")
gc()

ologit_model_lied42 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied4 + da_lied4 + ad_lied4 +
    politysenderminustarget + 
    lagtrade1 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied4, Hess = T,
  method = "logistic")
gc()

ologit_model_lied43 <- MASS::polr(
  as.factor(severity) ~ 
    aa_lied4 + da_lied4 + ad_lied4 +
    politysenderminustarget + 
    lagtrade3 +
    sender_hrs + target_hrs + 
    alliance + idealdistance + reviewer_revieweed + 
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_lied4, Hess = T,
  method = "logistic")
gc()

save(list = c(
  "ologit_model_poly0",  "ologit_model_poly1",  'ologit_model_poly2',  'ologit_model_poly3',
  'ologit_model_pol40',  'ologit_model_pol41',  'ologit_model_pol42',  'ologit_model_pol43',
  'ologit_model_brm0',   'ologit_model_brm1',   'ologit_model_brm2',   'ologit_model_brm3',
  'ologit_model_lied50', 'ologit_model_lied51', 'ologit_model_lied52', 'ologit_model_lied53',
  'ologit_model_lied40', 'ologit_model_lied41', 'ologit_model_lied42', 'ologit_model_lied43',
  'ologit_model_vdem0',  'ologit_model_vdem1',  'ologit_model_vdem2',  'ologit_model_vdem3'),
  file = "data/estimation/appendix_table_e1.RData")

#### Table E.1: Ordinal logit models on UPR recommendations shaming ------------
#### behaviors using different regime measurements
texreg::texreg(
  list(
    ologit_model_poly0, ologit_model_poly1, ologit_model_poly2, ologit_model_poly3,
    ologit_model_vdem0, ologit_model_vdem1, ologit_model_vdem2, ologit_model_vdem3,
    ologit_model_pol40, ologit_model_pol41, ologit_model_pol42, ologit_model_pol43,
    ologit_model_brm0,  ologit_model_brm1,  ologit_model_brm2,  ologit_model_brm3,
    ologit_model_lied40,ologit_model_lied41,ologit_model_lied42,ologit_model_lied43,
    ologit_model_lied50,ologit_model_lied51,ologit_model_lied52,ologit_model_lied53),
  omit.coef = "as.factor",
  file = "tables/appendix_table_e1.tex")
gc()

panels <- list(
  "EDI" = list(
    ologit_model_poly0, ologit_model_poly1, ologit_model_poly2, ologit_model_poly3),
  "RoW" = list(
    ologit_model_vdem0, ologit_model_vdem1, ologit_model_vdem2, ologit_model_vdem3),
  "POLITY5" = list(
    ologit_model_pol40, ologit_model_pol41, ologit_model_pol42, ologit_model_pol43),
  "BRM" = list(
    ologit_model_brm0,  ologit_model_brm1,  ologit_model_brm2,  ologit_model_brm3),
  "LIED < 4" = list(
    ologit_model_lied40,ologit_model_lied41,ologit_model_lied42,ologit_model_lied43),
  "LIED < 5" = list(
    ologit_model_lied50,ologit_model_lied51,ologit_model_lied52,ologit_model_lied53)
)
rows1 <- 
  tribble(~term,        
          ~ologit_model_poly0, ~ologit_model_poly1, ~ologit_model_poly2, ~ologit_model_poly3,
          ~ologit_model_vdem0, ~ologit_model_vdem1, ~ologit_model_vdem2, ~ologit_model_vdem3,
          ~ologit_model_pol40, ~ologit_model_pol41, ~ologit_model_pol42, ~ologit_model_pol43,
          ~ologit_model_brm0,  ~ologit_model_brm1,  ~ologit_model_brm2,  ~ologit_model_brm3,
          ~ologit_model_lied40,~ologit_model_lied41,~ologit_model_lied42,~ologit_model_lied43,
          ~ologit_model_lied50,~ologit_model_lied51,~ologit_model_lied52,~ologit_model_lied53,
          'Controls', 
          'No', 'Share of Combined GDP', 'L.Bilateral Trade$_\\text{t-1}$', 'L.Bilateral Trade$_\\text{t-3}$', 
          'No', 'Share of Combined GDP', 'L.Bilateral Trade$_\\text{t-1}$', 'L.Bilateral Trade$_\\text{t-3}$',
          'No', 'Share of Combined GDP', 'L.Bilateral Trade$_\\text{t-1}$', 'L.Bilateral Trade$_\\text{t-3}$',
          'No', 'Share of Combined GDP', 'L.Bilateral Trade$_\\text{t-1}$', 'L.Bilateral Trade$_\\text{t-3}$',
          'No', 'Share of Combined GDP', 'L.Bilateral Trade$_\\text{t-1}$', 'L.Bilateral Trade$_\\text{t-3}$',
          'No', 'Share of Combined GDP', 'L.Bilateral Trade$_\\text{t-1}$', 'L.Bilateral Trade$_\\text{t-3}$')

attr(rows1, 'position') <- c(7)
options("modelsummary_format_numeric_latex" = "plain")
modelsummary::modelsummary(
  estimate = "{estimate}{stars} ({std.error})",
  panels, stars = TRUE, shape = "rbind",
  coef_map = c("aa_poly" = "Autocracy Dyads",
               "ad_poly" = "Auto$_\\text{Reviewer}$-Demo$_\\text{SuR}$ Mixed Dyads",
               "da_poly" = "Demo$_\\text{Reviewer}$-Auto$_\\text{SuR}$ Mixed Dyads",
               "aa_vdem" = "Autocracy Dyads",
               "ad_vdem" = "Auto$_\\text{Reviewer}$-Demo$_\\text{SuR}$ Mixed Dyads",
               "da_vdem" = "Demo$_\\text{Reviewer}$-Auto$_\\text{SuR}$ Mixed Dyads",
               "aa_pol4" = "Autocracy Dyads",
               "ad_pol4" = "Auto$_\\text{Reviewer}$-Demo$_\\text{SuR}$ Mixed Dyads",
               "da_pol4" = "Demo$_\\text{Reviewer}$-Auto$_\\text{SuR}$ Mixed Dyads",
               "aa_brm" = "Autocracy Dyads",
               "ad_brm" = "Auto$_\\text{Reviewer}$-Demo$_\\text{SuR}$ Mixed Dyads",
               "da_brm" = "Demo$_\\text{Reviewer}$-Auto$_\\text{SuR}$ Mixed Dyads",
               "aa_lied4" = "Autocracy Dyads",
               "ad_lied4" = "Auto$_\\text{Reviewer}$-Demo$_\\text{SuR}$ Mixed Dyads",
               "da_lied4" = "Demo$_\\text{Reviewer}$-Auto$_\\text{SuR}$ Mixed Dyads",
               "aa_lied5" = "Autocracy Dyads",
               "ad_lied5" = "Auto$_\\text{Reviewer}$-Demo$_\\text{SuR}$ Mixed Dyads",
               "da_lied5" = "Demo$_\\text{Reviewer}$-Auto$_\\text{SuR}$ Mixed Dyads"),
  #add_rows = rows1,
  digits = 3, 
  coef_omit = "as.factor",
  gof_omit = 'DF|Deviance|R2|adjR2', 
  #output = "html"
  output = "tables/appendix_table_e1_stylized.tex"
)

### E.2 Alternative Explanation for Democratic Variants ------------------------
ERT::get_eps() -> ert
vdemdata::vdem -> vdem

vdem |> 
  dplyr::select(COWcode, year, v2xps_party, e_lexical_index) -> 
  vdem_party

ert |> 
  mutate(
    COWcode = 
      countrycode::countrycode(country_name, "country.name", "cown"),
    regime = if_else(v2x_regime < 2, 0L,
                     if_else(v2x_regime == 2L | 
                               v2x_regime == 3L, 1L, NA_integer_)),
    demoyear = 
      if_else(
        reg_trans == 1L & !is.na(year), as.integer(year),
        if_else(
          reg_trans < 1L & !is.na(year), 0L, NA_integer_)),
    demoyear2 = 
      if_else(
        reg_trans == 1L & !is.na(year), 1L,
        if_else(
          reg_trans < 1L & !is.na(year), 0L, NA_integer_))) ->
  ert_trans

ert_trans |> left_join(
  vdem_party |> dplyr::select(-v2xps_party) |> drop_na(),
  by = c("COWcode", "year")
) -> ert_vdem

ert_vdem |>
  mutate(regime = 
           if_else(e_lexical_index > 5, 1L,
                   if_else(e_lexical_index < 6, 0L, NA_integer_))) ->
  ert_vdem

ert_vdem |>
  group_by(COWcode, regime) |>
  mutate(duration = dplyr::row_number(),
         newdemo = if_else(duration > 30 & regime == 1, 1L,
                           if_else(duration <= 30 & regime == 1, 0L,
                                   NA_integer_))) |>
  ungroup() -> ert_vdem


sample |> mutate(to_cow = as.numeric(as.character(to_cow)),
                 from_cow = as.numeric(as.character(from_cow)),
                 year = as.numeric(as.character(year))) |> 
  left_join(ert_vdem |> 
              dplyr::select(COWcode, year, newdemo) |>
              drop_na(), 
            by = c("to_cow" = "COWcode", "year")) |>
  rename(target_newdemoid = newdemo) -> upr_party

upr_party |>
  left_join(ert_vdem |> 
              dplyr::select(COWcode, year, newdemo) |>
              drop_na(),
            by = c("from_cow" = "COWcode", "year")) |>
  rename(sender_newdemoid = newdemo) -> upr_party

upr_party |>
  left_join(vdem_party |> drop_na(), 
            by = c("to_cow" = "COWcode", "year")) |>
  rename(target_party = v2xps_party) -> upr_party

upr_party |>
  left_join(vdem_party |> drop_na(), 
            by = c("from_cow" = "COWcode", "year")) |>
  rename(sender_party = v2xps_party) -> upr_party

upr_party |>
  left_join(ert_trans |> drop_na(),
            by = c("to_cow" = "COWcode", "year")) |>
  rename(target_trans = reg_trans) -> upr_party_trans

upr_party_trans |>
  left_join(ert_trans |> drop_na(),
            by = c("from_cow" = "COWcode", "year")) |>
  rename(sender_trans = reg_trans) -> upr_party_trans

upr_party |>
  mutate(
    dyads_and = case_when(
      target_lexi_regime == "Autocracy" &
        sender_lexi_regime == "Autocracy" ~ "Autocratic Dyads",
      target_lexi_regime == "Autocracy" &
        sender_lexi_regime == "Democracy" &
        sender_newdemoid == 1L ~ "Auto(SUR)-New Democratic(REV) Dyads",
      target_lexi_regime == "Democracy" & 
        target_newdemoid == 1L &
        sender_lexi_regime == "Autocracy" ~ "New Democracy(SUR)-Auto(REV) Dyads",
      target_lexi_regime == "Democracy" & 
        target_newdemoid == 1L &
        sender_lexi_regime == "Democracy" &
        sender_newdemoid == 1L ~ "New Democratic Dyads",
      T ~ NA_character_
    ),
    dyads_aod = case_when(
      target_lexi_regime == "Autocracy" &
        sender_lexi_regime == "Autocracy" ~ "Autocratic Dyads",
      target_lexi_regime == "Autocracy" &
        sender_lexi_regime == "Democracy" &
        sender_newdemoid == 0L ~ "Auto(SUR)-Old Democratic(REV) Dyads",
      target_lexi_regime == "Democracy" & 
        target_newdemoid == 0L &
        sender_lexi_regime == "Autocracy" ~ "Old Democracy(SUR)-Auto(REV) Dyads",
      target_lexi_regime == "Democracy" & 
        target_newdemoid == 0L &
        sender_lexi_regime == "Democracy" &
        sender_newdemoid == 0L ~ "Old Democratic Dyads",
      T ~ NA_character_
    )) -> upr_party_dyads

sample |> group_by(rowid) |> count() |> dplyr::filter(n > 1) -> dup

sample |> dplyr::filter(!rowid %in% dup$rowid) |> 
  left_join(
    upr_party_dyads |> dplyr::select(rowid, dyads_and, dyads_aod) |> 
      dplyr::filter(!rowid %in% dup$rowid),
    by = "rowid"
  ) -> newold

newold |> dplyr::select(
  severity, 
  dyads_and,
  share_combined,
  politysenderminustarget,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_regime_and

sample_regime_and <- as.data.frame(sample_regime_and)

sample_regime_and |>
  mutate(
    ndnd = if_else(dyads_and == "New Democratic Dyads", 1L, 0L),
    aa = if_else(dyads_and == "Autocratic Dyads", 1L, 0L),
    and = if_else(dyads_and == "Auto(SUR)-New Democratic(REV) Dyads", 1L, 0L),
    nda = if_else(dyads_and == "New Democracy(SUR)-Auto(REV) Dyads", 1L, 0L),
    dyads_and = factor(dyads_and, 
                       levels = c(
                         "Autocratic Dyads", 
                         "New Democratic Dyads",
                         "Auto(SUR)-New Democratic(REV) Dyads",
                         "NNew Democracy(SUR)-Auto(REV) Dyads"))
  ) -> sample_regime_and

newold |> dplyr::select(
  severity, 
  dyads_aod,
  share_combined,
  politysenderminustarget,
  sender_hrs, target_hrs,
  lagtrade1, lagtrade3,
  alliance, idealdistance, reviewer_revieweed, 
  to_mem, from_mem, regionalism, 
  migrantscat, socioecon, vulnerable, 
  womenchild, physint, justicecat, political, 
  to_cow, from_cow, 
  year) |> drop_na(severity) -> sample_regime_aod

sample_regime_aod <- as.data.frame(sample_regime_aod)

sample_regime_aod |>
  mutate(
    odod = if_else(dyads_aod == "Old Democratic Dyads", 1L, 0L),
    aa = if_else(dyads_aod == "Autocratic Dyads", 1L, 0L),
    aod = if_else(dyads_aod == "Auto(SUR)-Old Democratic(REV) Dyads", 1L, 0L),
    oda = if_else(dyads_aod == "Old Democracy(SUR)-Auto(REV) Dyads", 1L, 0L),
    dyads_aod = factor(dyads_aod, 
                       levels = c("Autocratic Dyads", 
                                  "Auto(SUR)-Old Democratic(REV) Dyads",
                                  "Old Democracy(SUR)-Auto(REV) Dyads",
                                  "Old Democratic Dyads"))
  ) -> sample_regime_aod

sample_regime_and$dyads_and <- factor(
  sample_regime_and$dyads_and,
  levels = c("Autocratic Dyads", "Autocracy(SUR)-New Democracy(REV) Dyads",
             "New Democracy(SUR)-Autocracy(REV) Dyads", "New Democracy Dyads")
)
sample_regime_aod$dyads_aod <- factor(
  sample_regime_aod$dyads_aod,
  levels = c("Autocratic Dyads", "Auto(SUR)-Old Democratic(REV) Dyads",
             "Old Democracy(SUR)-Auto(REV) Dyads", "Old Democratic Dyads")
)
gc()

sample_regime_and |> mutate(severity = factor(severity, levels = c(1, 2, 3))) ->
  sample_regime_and
sample_regime_aod |> mutate(severity = factor(severity, levels = c(1, 2, 3))) ->
  sample_regime_aod

ologit_model_and1 <- MASS::polr(
  as.factor(severity) ~ 
    ndnd + and + nda +
    politysenderminustarget +
    share_combined +
    #log(abs(gdpgap+1)) +
    alliance + idealdistance + reviewer_revieweed + 
    sender_hrs + target_hrs +
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_regime_and,
  method = "logistic", Hess = T)
gc()

ologit_model_and2 <- MASS::polr(
  as.factor(severity) ~ 
    ndnd + and + nda +
    politysenderminustarget + 
    lagtrade1 +
    alliance + idealdistance + reviewer_revieweed + 
    sender_hrs + target_hrs +
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_regime_and,
  Hess = T,
  method = "logistic")
gc()


ologit_model_aod1 <- MASS::polr(
  as.factor(severity) ~ 
    aod + oda  + odod + 
    politysenderminustarget + 
    share_combined +
    #log(abs(gdpgap+1)) +
    alliance + idealdistance + reviewer_revieweed + 
    sender_hrs + target_hrs +
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_regime_aod,
  Hess = T,
  method = "logistic")
gc()

ologit_model_aod2 <- MASS::polr(
  as.factor(severity) ~ 
    aod + oda  + odod + 
    politysenderminustarget + 
    lagtrade1 +
    alliance + idealdistance + reviewer_revieweed + 
    sender_hrs + target_hrs +
    to_mem + from_mem + regionalism +
    migrantscat + socioecon + vulnerable + 
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_regime_aod,
  Hess = T,
  method = "logistic")
gc()

#### Table E.2: Ordinal logit models on UPR recommendations shaming ------------
#### behaviors by different types of democracy
texreg::texreg(list(
  ologit_model_and1,ologit_model_and2,
  ologit_model_aod1, ologit_model_aod2),
  omit.coef = "as.factor",
  single.row = T,
  custom.coef.map=list(
    "ndnd" = "New Democracy Dyads", 
    "and"  =  "Auto$_{\\text{SuR}}$-New Demo$_{\\text{Reviewer}}$ Dyads",
    "nda"  =  "New Demo$_{\\text{SuR}}$-Auto$_{\\text{Reviewer}}$ Dyads", 
    "politysenderminustarget" = "EDI difference$_{\\text{Reveiwer - State under review}}$",
    "share_combined" = "Share of Combined GDP",
    "alliance" = "Alliance",
    "idealdistance" = "Geopolitical Affinity",
    "reviewer_revieweed" = "Reviewed Reviewer",
    "sender_hrs" = "Human Rights Score$_{\\text{Reviewer}}$",
    "target_hrs" = "Human Rights Score$_{\\text{SuR}}$",
    "to_mem" = "HRC Member$_{\\text{SuR}}$", 
    "from_mem" = "HRC Member$_{\\text{Reviewer}}$",
    "regionalism" = "Same Region",
    "migrantscat" = "Migration",
    "socioecon" = "Socio-Economic Rights", "vulnerable" = "Vulnerable Populations",
    "womenchild" = "Women, Children \\& Trafficking", "physint" = "Physical Integrity Rights",
    "justicecat" = "Justice", "political" = "Speech \\& Political Participation",
    "lagtrade1" = "L.ln(trade$_{t-1}$)",
    "oda" = "Old Demo$_{\\text{SuR}}$-Auto$_{\\text{Reviewer}}$ Dyads", 
    "aod" = "Auto$_{\\text{SuR}}$-Old Demo$_{\\text{Reviewer}}$ Dyads", "odod" = "Old Demo Dyads",
    "(Cut 1)", "(Cut 2)"),
  reorder.coef = c(1, 2, 3, 22, 23, 24, 5, 21, 4, 6, 7, 8, 9, 11, 
                   10, 12, 13, 14, 15, 16, 17, 18, 19, 20),
  include.rsquared=FALSE, 
  custom.gof.rows=list(`Country-fixed` = c("Yes", "Yes", "Yes", "Yes"),
                       `Year-fixed` = c("Yes", "Yes", "Yes", "Yes")),
  stars=c(0.001, 0.01, 0.05), 
  #custom.note="%stars. Year effects not significant.",
  label = "appendix-tbl-e2",
  caption=
    "Ordinal logit models on UPR recommendations shaming behaviors",
  file = "tables/appendix_table_e2.tex")

### First difference for new democracy dyads
#### Predicted Probabilities
##### Here, let's use the same profile as first difference estimation
##### Matrix to store results
set.seed(1234)
simb_nd <- mvrnorm(n = 4000, 
                   mu = c(coef(ologit_model_and1), ologit_model_and1$zeta),
                   Sigma = vcov(ologit_model_and1))

simb_od <- mvrnorm(n = 4000, 
                   mu = c(coef(ologit_model_aod1), ologit_model_aod1$zeta),
                   Sigma = vcov(ologit_model_aod1))


cut_nd <- cbind(-Inf, simb_nd[, 289:290], Inf)
cut_od <- cbind(-Inf, simb_od[, 275:276], Inf)
#### Create the baseline profile + Calculate the Predicted Probs
x_dd_nd <- c(0, 0, 0, 
             mean(sample_regime_and$politysenderminustarget, na.rm = T),
             mean(sample_regime_and$share_combined, na.rm = T),
             0, 
             mean(sample_regime_and$idealdistance, na.rm = T), 
             0, 
             mean(sample_regime_and$sender_hrs, na.rm = T), 
             mean(sample_regime_and$target_hrs, na.rm = T), 
             0, 0, 0, 
             0, 0, 0, 0, 0, 0, 0, rep(0, 270))
xddb_nd <- simb_nd[, 1:290] %*% x_dd_nd

pr0_nd <- matrix(NA, nrow = 4000, ncol = 3)
for (i in 1:3){ # for each value of the dependent variable...
  pr0_nd[, i] <- plogis(cut_nd[, i + 1] - xddb_nd) - plogis(cut_nd[, i] - xddb_nd)
}
#### Matrix to store the results
out.mat_nd <- matrix(NA, nrow = 3 * 3, ncol = 3)
rownames(out.mat_nd) <- paste0(rep(c("New Democracy Dyads",
                                     "Auto-New Demo Mixed Dyads",
                                     "New Demo-Auto Mixed Dyads"), 
                                   each = 3), "-",
                               rep(1:3, times = 3))
colnames(out.mat_nd) <- c("Mean", "Lower", "Upper")

#### Then compute the first differences in each covariates
#### DD -> AA
x_aa_nd <- x_dd_nd
x_aa_nd[1] <- 1

x_aa_ndb <- simb_nd[, 1:290] %*% x_aa_nd
for (i in 1:3){ # for each value of the dependent variable...
  pr1_nd <- plogis(cut_nd[, i + 1] - x_aa_ndb) - plogis(cut_nd[, i] - x_aa_ndb)
  fdiff <- pr1_nd - pr0_nd[, i]
  out.mat_nd[i, 1] <- mean(fdiff)
  out.mat_nd[i, 2:3] <- quantile(fdiff, probs = c(0.025, 0.975))
}

#### DD -> DA
x_da_nd <- x_dd_nd
x_da_nd[2] <- 1
xdab <- simb_nd[, 1:290] %*% x_da_nd
for (i in 1:3){ # for each value of the dependent variable...
  pr1_nd <- plogis(cut_nd[, i + 1] - xdab) - plogis(cut_nd[, i] - xdab)
  fdiff <- pr1_nd - pr0_nd[, i]
  out.mat_nd[i + 3, 1] <- mean(fdiff)
  out.mat_nd[i + 3, 2:3] <- quantile(fdiff, probs = c(0.025, 0.975))
}

#### DD -> AD
x_ad_nd <- x_dd_nd
x_ad_nd[3] <- 1
xadb <- simb_nd[, 1:290] %*% x_ad_nd
for (i in 1:3){ # for each value of the dependent variable...
  pr1_nd <- plogis(cut_nd[, i + 1] - xadb) - plogis(cut_nd[, i] - xadb)
  fdiff <- pr1_nd - pr0_nd[, i]
  out.mat_nd[i + 6, 1] <- mean(fdiff)
  out.mat_nd[i + 6, 2:3] <- quantile(fdiff, probs = c(0.025, 0.975))
}
out.mat_nd # print the simulation results

# Visualize the results
library(wesanderson)
out.mat_nd <- as.data.frame(out.mat_nd)
out.mat_nd |> mutate(
  sig = if_else(Lower < 0 & Upper > 0, 
                "Insignificant", "Significant")
) -> out.mat_nd
out.mat_nd$sig <- factor(out.mat_nd$sig, 
                         levels = c("Significant",
                                    "Insignificant"))
out.mat_nd$regimedyads <- rep(1:3, times = 3)
g1_nd <- ggplot(out.mat_nd[1:3,]) +
  geom_pointrange(aes(x = regimedyads, y = Mean, 
                      ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F) +
  labs(subtitle = "Autocratic Dyads\nto New Democracy Dyads",
       y = "First Differences in Predicted Probabilities",
       x = NULL) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.35, 0.3)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  # ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean,
  #                              label = paste0(round(Mean, 4)*100, "%"),
  #                              color = sig),
  #                          size = 2.5, angle = 90,
  #                         # nudge_x = -0.1,
  #                          xlim = c(1, 3), #ylim = c(-Inf, Inf),
  #                          show.legend = F) +
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[1],
                                wesanderson::wes_palette("Royal1")[2])) + 
  theme(legend.title = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"),
        legend.position="none")

g1_nd
g2_nd <- ggplot(out.mat_nd[4:6,]) +
  geom_pointrange(aes(x = regimedyads, y = Mean, 
                      ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F) +
  labs(subtitle = "Autocracy Dyads\nto New Demo(SUR)-Auto(REV) Mixed Dyads",
       y = NULL, x = NULL) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.35, 0.3)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  # ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean,
  #                              label = paste0(round(Mean, 4)*100, "%"),
  #                              color = sig),
  #                         # nudge_x = -0.1,
  #                          size = 2.5, angle = 90,
  #                          xlim = c(1, 3), #ylim = c(-Inf, Inf),
  #                          show.legend = F) +
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[2],
                                wesanderson::wes_palette("Royal1")[1])) + 
  theme(legend.title = element_blank(),
        axis.text.y.left = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"),
        legend.position="none")
g2_nd
g3_nd <- ggplot(out.mat_nd[7:9,]) +
  geom_pointrange(aes(x = regimedyads, 
                      y = Mean, ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F) +
  labs(subtitle = "Autocracy Dyads\nto Auto(SUR)-New Demo Mixed Dyads",
       x = NULL, y = NULL) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.35, 0.3)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  # ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean,
  #                              label = paste0(round(Mean, 4)*100, "%"),
  #                              color = sig),
  #                          size = 2.5, angle = 90,
  #                        #  nudge_x = -0.1,
  #                          xlim = c(1, 3), #ylim = c(-Inf, Inf),
  #                          show.legend = F) +
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[1],
                                wesanderson::wes_palette("Royal1")[2])) + 
  theme(legend.title = element_blank(),
        axis.text.y.left = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"),
        legend.position="none")

g3_nd

g1_nd + g2_nd + g3_nd + patchwork::plot_layout(ncol = 3)  -> Appendix_ologit_first_nd

#### Figure E.1: The first difference change in predicted probability by regime dyads: -----
####             New democracy and autocracy
Appendix_ologit_first_nd
ggsave("figures/appendix_figure_e1.pdf", 
       height = 4, width = 10, dpi = 900)

### First difference for old democracy dyads
x_dd_od <- c(0, 0, 0, 
             mean(sample_regime_and$politysenderminustarget, na.rm = T),
             mean(sample_regime_and$share_combined, na.rm = T),
             0, 
             mean(sample_regime_and$idealdistance, na.rm = T), 
             0, 
             mean(sample_regime_and$sender_hrs, na.rm = T),
             mean(sample_regime_and$target_hrs, na.rm = T), 
             0, 0, 0, 
             0, 0, 0, 0, 0, 0, 0, rep(0, 256))

xddb_od <- simb_od[, 1:276] %*% x_dd_od

pr0_od <- matrix(NA, nrow = 4000, ncol = 3)
for (i in 1:3){ # for each value of the dependent variable...
  pr0_od[, i] <- plogis(cut_od[, i + 1] - xddb_od) - plogis(cut_od[, i] - xddb_od)
}
#### Matrix to store the results
out.mat_od <- matrix(NA, nrow = 3 * 3, ncol = 3)
rownames(out.mat_od) <- paste0(rep(c("Auto(SUR)-Old Demo(REV) Mixed Dyads",
                                     "Old Demo(SUR)-Auto Mixed(REV) Dyads",
                                     "Old Democracy Dyads"), 
                                   each = 3), "-",
                               rep(1:3, times = 3))
colnames(out.mat_od) <- c("Mean", "Lower", "Upper")

#### Then compute the first differences in each covariates
#### DD -> AA
x_aa_od <- x_dd_od
x_aa_od[1] <- 1

x_aa_odb <- simb_od[, 1:276] %*% x_aa_od
for (i in 1:3){ # for each value of the dependent variable...
  pr1_od <- plogis(cut_od[, i + 1] - x_aa_odb) - plogis(cut_od[, i] - x_aa_odb)
  fdiff <- pr1_od - pr0_od[, i]
  out.mat_od[i, 1] <- mean(fdiff)
  out.mat_od[i, 2:3] <- quantile(fdiff, probs = c(0.025, 0.975))
}

#### DD -> DA
x_da_od <- x_dd_od
x_da_od[2] <- 1
xdab <- simb_od[, 1:276] %*% x_da_od
for (i in 1:3){ # for each value of the dependent variable...
  pr1_od <- plogis(cut_od[, i + 1] - xdab) - plogis(cut_od[, i] - xdab)
  fdiff <- pr1_od - pr0_od[, i]
  out.mat_od[i + 3, 1] <- mean(fdiff)
  out.mat_od[i + 3, 2:3] <- quantile(fdiff, probs = c(0.025, 0.975))
}

#### DD -> AD
x_ad_od <- x_dd_od
x_ad_od[3] <- 1
xadb <- simb_od[, 1:276] %*% x_ad_od
for (i in 1:3){ # for each value of the dependent variable...
  pr1_od <- plogis(cut_od[, i + 1] - xadb) - plogis(cut_od[, i] - xadb)
  fdiff <- pr1_od - pr0_od[, i]
  out.mat_od[i + 6, 1] <- mean(fdiff)
  out.mat_od[i + 6, 2:3] <- quantile(fdiff, probs = c(0.025, 0.975))
}
out.mat_od # print the simulation results

# Visualize the results
library(wesanderson)
out.mat_od <- as.data.frame(out.mat_od)
out.mat_od |> mutate(
  sig = if_else(Lower < 0 & Upper > 0, 
                "Insignificant", "Significant")
) -> out.mat_od
out.mat_od$sig <- factor(out.mat_od$sig, 
                         levels = c("Significant",
                                    "Insignificant"))
out.mat_od$regimedyads <- rep(1:3, times = 3)
g1_od <- ggplot(out.mat_od[1:3,]) +
  geom_pointrange(aes(x = regimedyads, y = Mean, 
                      ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F) +
  labs(subtitle = "Autocracy Dyads\nto Auto(SUR)-Old Demo(REV) Mixed Dyads",
       y = "First Differences in Predicted Probabilities",
       x = NULL) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.35, 0.35)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  # ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean,
  #                              label = paste0(round(Mean, 4)*100, "%"),
  #                              color = sig),
  #                          size = 2.5, angle = 90,
  #                         # nudge_x = -0.1,
  #                          xlim = c(1, 3), #ylim = c(-Inf, Inf),
  #                          show.legend = F) +
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[2],
                                wesanderson::wes_palette("Royal1")[1])) + 
  theme(legend.title = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"),
        legend.position="none")

g1_od
g2_od <- ggplot(out.mat_od[4:6,]) +
  geom_pointrange(aes(x = regimedyads, y = Mean, 
                      ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F) +
  labs(subtitle = "Autocracy Dyads\nto Old Demo(SUR)-Auto(REV) Mixed Dyads",
       y = NULL, x = NULL) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.35, 0.35)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  # ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean,
  #                              label = paste0(round(Mean, 4)*100, "%"),
  #                              color = sig),
  #                         # nudge_x = -0.1,
  #                          size = 2.5, angle = 90,
  #                          xlim = c(1, 3), #ylim = c(-Inf, Inf),
  #                          show.legend = F) +
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[2],
                                wesanderson::wes_palette("Royal1")[1])) + 
  theme(legend.title = element_blank(),
        axis.text.y.left = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"),
        legend.position="none")
g2_od
g3_od <- ggplot(out.mat_od[7:9,]) +
  geom_pointrange(aes(x = regimedyads, 
                      y = Mean, ymin = Lower, ymax = Upper,
                      color = sig), show.legend = F) +
  labs(subtitle = "Autocracy Dyads\nto Old Democracy Dyads",
       x = NULL, y = NULL) + 
  geom_hline(aes(yintercept = 0), linetype = "dashed") +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(-0.35, 0.35)) +
  scale_x_continuous(breaks = c(1:3),
                     labels = c("Praising",
                                "Neutral",
                                "Shaming")) + 
  # ggrepel::geom_text_repel(aes(x = regimedyads, y = Mean,
  #                              label = paste0(round(Mean, 4)*100, "%"),
  #                              color = sig),
  #                          size = 2.5, angle = 90,
  #                        #  nudge_x = -0.1,
  #                          xlim = c(1, 3), #ylim = c(-Inf, Inf),
  #                          show.legend = F) +
  scale_color_manual(values = c(wesanderson::wes_palette("Royal1")[1],
                                wesanderson::wes_palette("Royal1")[2])) + 
  theme(legend.title = element_blank(),
        axis.text.y.left = element_blank(),
        plot.margin = unit(c(0,1,0,0), "cm"),
        legend.position="none")

g3_od

g1_od + g2_od + g3_od + patchwork::plot_layout(ncol = 3)  ->
  Appendix_ologit_first_od


#### Figure E.2: The first difference change in predicted probability by regime dyads: -----
####             Old democracy and autocracy
Appendix_ologit_first_od
ggsave("figures/appendix_figure_e2.pdf", height = 4, width = 10, dpi = 900)

### E.3 Alternative Explanation for Autocratic Variants ------------------------
#### No China Sample ----------------------------------------------------

ologit_model1_ch <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 710), Hess = T,
  method = "logistic")
gc()

ologit_model2_ch <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    #ln2gdpgap +
    #ln_combined_gdppc +
    share_combined +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 710), Hess = T,
  method = "logistic")
gc()

ologit_model3_ch <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade1 +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 710), Hess = T,
  method = "logistic")
gc()

ologit_model4_ch <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade3 +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 710), Hess = T,
  method = "logistic")
gc()

#### No Russia Sample ----------------------------------------------------
ologit_model1_rs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 365), Hess = T,
  method = "logistic")

ologit_model2_rs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    #ln2gdpgap +
    #ln_combined_gdppc +
    share_combined +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 365), Hess = T,
  method = "logistic")
gc()

ologit_model3_rs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade1 +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 365), Hess = T,
  method = "logistic")
gc()

ologit_model4_rs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade3 +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% 365), Hess = T,
  method = "logistic")
gc()

#### No China or Russia Sample ----------------------------------------------------
ologit_model1_chrs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% c(365, 710)), Hess = T,
  method = "logistic")

ologit_model2_chrs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    #ln2gdpgap +
    #ln_combined_gdppc +
    share_combined +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% c(365, 710)), Hess = T,
  method = "logistic")
gc()

ologit_model3_chrs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade1 +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% c(365, 710)), Hess = T,
  method = "logistic")
gc()

ologit_model4_chrs <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    lagtrade3 +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    as.factor(to_cow) + as.factor(from_cow) +
    as.factor(year),
  data = sample_df |> dplyr::filter(!to_cow %in% c(710, 365)), Hess = T,
  method = "logistic")
gc()

save(ologit_model1_ch, ologit_model2_ch, ologit_model3_ch, ologit_model4_ch,
     ologit_model1_rs, ologit_model2_rs, ologit_model3_rs, ologit_model4_rs,
     ologit_model1_chrs, ologit_model2_chrs, ologit_model3_chrs, ologit_model4_chrs,
     file = "data/estimation/ologit_china_russia.rdata")

#### Table E.3: Ordinal logit models on UPR recommendations shaming behaviors ----
####            without Russia, China, or Both
texreg::texreg(list( 
  ologit_model1_ch, ologit_model2_ch, ologit_model3_ch, ologit_model4_ch,
  ologit_model1_rs, ologit_model2_rs, ologit_model3_rs, ologit_model4_rs,
  ologit_model1_chrs, ologit_model2_chrs, ologit_model3_chrs, ologit_model4_chrs),
  omit.coef = "as.factor",
  custom.model.names=c(
    "Model 1", "Model 2", "Model 3", "Model 4",
    "Model 5", "Model 6", "Model 7", "Model 8",
    "Model 9", "Model 10", "Model 11", "Model 12"),
  custom.coef.map = list(
    "aa_lexi" = "Autocracy Dyads",
    "da_lexi" = "Demo$_{\\mathrm{Reviewer}}$-Auto$_{\\mathrm{SuR}}$ Mixed Dyads", 
    "ad_lexi" = "Auto$_{\\mathrm{Reviewer}}$-Demo(SUR)",
    "1|2" = "(Cut 1)",
    "2|3" = "(Cut 2)",
    "politysenderminustarget" = "Joint Democracy", 
    "share_combined" = "Share of Combined GDP",
    "alliance" = "Alliance", 
    "idealdistance" =  "Geopolitical Affinity",
    "reviewer_revieweed" = "Reviewed Reviewer",
    "to_mem" = "HRC Member$_{\\mathrm{SuR}}$", 
    "from_mem" = "HRC Member$_{\\mathrm{Reviewer}}$",
    "regionalism"= "Same Region",
    "sender_hrs" = "Human Rights Score$_{\\mathrm{Reviewer}}$",
    "target_hrs" = "Human Rights Score$_{\\mathrm{SuR}}$",
    "migrantscat"= "Migration", 
    "socioecon" = "Socio-Economic Rights", 
    "vulnerable" = "Vulnerable Populations",
    "womenchild" =   "Women, Children \\& Trafficking", 
    "physint" = "Physical Integrity Rights",
    "justicecat" =   "Justice", 
    "political" = "Speech \\& Political Participation",
    "lagtrade1" =   "L.ln(trade$_{t−1}$)",
    "lagtrade3" =   "L.ln(trade$_{t−3}$)"),
  reorder.coef = c(1, 2, 3, 7, 23, 24, 6, 8, 9, 10, 11, 12, 13, 14,
                   15, 16, 17, 18, 19, 20, 21, 22, 4, 5),
  single.row = FALSE,
  include.rsquared=FALSE, 
  custom.gof.rows=list(`Country-fixed` = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                       `Year-fixed` = c("Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes")),
  stars=c(0.001, 0.01, 0.05), 
  #custom.note="%stars. Year effects not significant.",
  caption=
    "\\label{tab1} Ordered logit models on UPR recommendations shaming behaviors",
  file = "tables/appendix_table_e3.tex")

#### Figure E.3: Differences in simulated coefficients between ------------------
####             main results and sub-sample analyses for the democratic 
####             Reviewer and autocratic States under review mixed dyads on
####             recommendation specificity
ologit_previous <- MASS::polr(
  as.factor(severity) ~ 
    aa_lexi + da_lexi + ad_lexi +
    politysenderminustarget + 
    share_combined +
    alliance + idealdistance + reviewer_revieweed +
    to_mem + from_mem + regionalism +
    sender_hrs + target_hrs +
    migrantscat + socioecon + vulnerable +
    womenchild + physint + justicecat + political +
    to_cow + from_cow +
    year,
  data = sample,
  Hess = T, 
  method = "logistic")
gc()

### Main results
set.seed(1234)
simb_pre <- 
  mvrnorm(n = 4000, 
          mu = c(coef(ologit_previous), ologit_previous$zeta),
          Sigma = vcov(ologit_previous))

simb_pre[, 1:3] |> as_tibble() |> 
  tidyr::gather(`Regime dyads`, values) |> 
  group_by(`Regime dyads`) |> 
  summarize(
    mean = mean(values, na.rm = T),
    ll90 = quantile(values, probs = 0.05),
    ll95 = quantile(values, probs = 0.025),
    ul90 = quantile(values, probs = 0.95),
    ul95 = quantile(values, probs = 0.975)
  ) |> mutate(Model = "Previous results") -> pre

tibble(
  `Regime dyads` = c("aa_lexi", "da_lexi", "ad_lexi"),
  mean = c(-0.74, -0.16, 0.05),
  ll90 = c(-0.74 - 1.645*0.08, -0.16 - 1.645*0.07, 0.05 - 1.645*0.04),
  ll95 = c(-0.74 - 1.96*0.08, -0.16 - 1.96*0.07, 0.05 - 1.96*0.04),
  ul90 = c(-0.74 + 1.645*0.08, -0.16 + 1.645*0.07, 0.05 + 1.645*0.04),
  ul95 = c(-0.74 + 1.96*0.08, -0.16 + 1.96*0.07, 0.05 + 1.96*0.04),
  Model = c("Previous results", "Previous results", "Previous results")
) -> pre

simb <-   mvrnorm(n = 4000, 
                  mu = c(coef(ologit_model2), ologit_model2$zeta),
                  Sigma = vcov(ologit_model2))

simb[, 1:3] |> as_tibble() |> 
  tidyr::gather(`Regime dyads`, values) |> 
  group_by(`Regime dyads`) |> 
  summarize(
    mean = mean(values, na.rm = T),
    ll90 = quantile(values, probs = 0.05),
    ll95 = quantile(values, probs = 0.025),
    ul90 = quantile(values, probs = 0.95),
    ul95 = quantile(values, probs = 0.975)
  ) |> mutate(Model = "Main results") -> main

#### Russia/China/Both
set.seed(1234)

simb_norus <- 
  mvrnorm(n = 4000, 
          mu = c(coef(ologit_model1_rs), ologit_model1_rs$zeta),
          Sigma = vcov(ologit_model1_rs))

simb_norus[, 1:3] |> as_tibble() |> 
  tidyr::gather(`Regime dyads`, values) |> 
  group_by(`Regime dyads`) |> 
  summarize(
    mean = mean(values, na.rm = T),
    ll90 = quantile(values, probs = 0.05),
    ll95 = quantile(values, probs = 0.025),
    ul90 = quantile(values, probs = 0.95),
    ul95 = quantile(values, probs = 0.975)
  ) |> mutate(Model = "W/o Russia") -> norussia

set.seed(1234)
simb_nochina <- 
  mvrnorm(n = 4000, 
          mu = c(coef(ologit_model1_ch), ologit_model1_ch$zeta),
          Sigma = vcov(ologit_model1_ch))

simb_nochina[, 1:3] |> as_tibble() |> 
  tidyr::gather(`Regime dyads`, values) |> 
  group_by(`Regime dyads`) |> 
  summarize(
    mean = mean(values, na.rm = T),
    ll90 = quantile(values, probs = 0.05),
    ll95 = quantile(values, probs = 0.025),
    ul90 = quantile(values, probs = 0.95),
    ul95 = quantile(values, probs = 0.975)
  ) |> mutate(Model = "W/o China") -> nochina

set.seed(1234)
simb_noboth <- 
  mvrnorm(n = 4000, 
          mu = c(coef(ologit_model1_chrs), ologit_model1_chrs$zeta),
          Sigma = vcov(ologit_model1_chrs))

simb_noboth[, 1:3] |> as_tibble() |> 
  tidyr::gather(`Regime dyads`, values) |> 
  group_by(`Regime dyads`) |> 
  summarize(
    mean = mean(values, na.rm = T),
    ll90 = quantile(values, probs = 0.05),
    ll95 = quantile(values, probs = 0.025),
    ul90 = quantile(values, probs = 0.95),
    ul95 = quantile(values, probs = 0.975)
  ) |> mutate(Model = "W/o Both") -> noboth

tibble(
  Model = c("W/o Both", "W/o Both", "W/o Both",
            "W/o China", "W/o China", "W/o China",
            "W/o Russia", "W/o Russia", "W/o Russia"),
  Dyad = c("Autocratic Dyads: H1", 
           "Demo (Reviewer)-Auto (SuR) Mixed Dyads: H3",
           "Auto (Reviewer)-Demo (SuR) Mixed Dyads: H2",
           "Autocratic Dyads: H1", 
           "Demo (Reviewer)-Auto (SuR) Mixed Dyads: H3",
           "Auto (Reviewer)-Demo (SuR) Mixed Dyads: H2",
           "Autocratic Dyads: H1", 
           "Demo (Reviewer)-Auto (SuR) Mixed Dyads: H3",
           "Auto (Reviewer)-Demo (SuR) Mixed Dyads: H2"),
  mean = c(as.numeric(main[1, 2] - noboth[1, 2]),
           as.numeric(main[2, 2] - noboth[2, 2]),
           as.numeric(main[3, 2] - noboth[3, 2]),
           as.numeric(main[1, 2] - nochina[1, 2]),
           as.numeric(main[2, 2] - nochina[2, 2]),
           as.numeric(main[3, 2] - nochina[3, 2]),
           as.numeric(main[1, 2] - norussia[1, 2]),
           as.numeric(main[2, 2] - norussia[2, 2]),
           as.numeric(main[3, 2] - norussia[3, 2])),
  ll = c(as.numeric(main[1, 4] - noboth[1, 4]),
         as.numeric(main[2, 4] - noboth[2, 4]),
         as.numeric(main[3, 4] - noboth[3, 4]),
         as.numeric(main[1, 4] - nochina[1, 4]),
         as.numeric(main[2, 4] - nochina[2, 4]),
         as.numeric(main[3, 4] - nochina[3, 4]),
         as.numeric(main[1, 4] - norussia[1, 4]),
         as.numeric(main[2, 4] - norussia[2, 4]),
         as.numeric(main[3, 4] - norussia[3, 4])),
  ul = c(as.numeric(main[1, 5] - noboth[1, 5]),
         as.numeric(main[2, 5] - noboth[2, 5]),
         as.numeric(main[3, 5] - noboth[3, 5]),
         as.numeric(main[1, 5] - nochina[1, 5]),
         as.numeric(main[2, 5] - nochina[2, 5]),
         as.numeric(main[3, 5] - nochina[3, 5]),
         as.numeric(main[1, 5] - norussia[1, 5]),
         as.numeric(main[2, 5] - norussia[2, 5]),
         as.numeric(main[3, 5] - norussia[3, 5])))

tibble(
  Model = rep(c(rep("W/o China", 4000), rep("W/o Russia", 4000), 
                rep("W/o Both", 4000)), 3),
  Diff = c(simb_nochina[, "aa_lexi"] - simb[, "aa_lexi"],
           simb_norus[, "aa_lexi"] - simb[, "aa_lexi"],
           simb_noboth[, "aa_lexi"] - simb[, "aa_lexi"],
           simb_nochina[, "da_lexi"] - simb[, "da_lexi"],
           simb_norus[, "da_lexi"] - simb[, "da_lexi"],
           simb_noboth[, "da_lexi"] - simb[, "da_lexi"],
           simb_nochina[, "ad_lexi"] - simb[, "ad_lexi"],
           simb_norus[, "ad_lexi"] - simb[, "ad_lexi"],
           simb_noboth[, "ad_lexi"] - simb[, "ad_lexi"]),
  Dyads = c(rep("Autocracy Dyads: H1", 12000),
            rep("Demo (Reviewer)-Auto (SuR)\nMixed Dyads: H3", 12000),
            rep("Auto (Reviewer)-Demo (SuR)\nMixed Dyads: H2", 12000))) |> 
  mutate(
    Dyads = 
      factor(Dyads,
             levels = c("Autocracy Dyads: H1",
                        "Auto (Reviewer)-Demo (SuR)\nMixed Dyads: H2",
                        "Demo (Reviewer)-Auto (SuR)\nMixed Dyads: H3"))) |> 
  ggplot(aes(Diff)) + 
  geom_density(color = "black", fill = "white") + 
  geom_histogram(aes(y = ..density..), color = "black", fill = "white") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "darkred") +
  labs(x = "\nDifferences in estimated coefficients",
       y = "Density\n", 
       caption = str_wrap("Note: Each panel represents the distribution of differences in simulated coefficients for different dyads: Autocracy Dyads (H1), Demo (Reviewer) - Auto (SuR) Mixed Dyads (H2), and Auto (Reviewer) - Demo (SuR) Mixed Dyads (H3). The vertical dashed red line at 0 helps to visualize the point where there is no difference between the main and sub-sample models.",
                          115, exdent = 10)) +
  facet_wrap(~Model + Dyads) +
  theme(
    legend.title = element_blank(),
    legend.text = element_text(size = rel(1.3)),
    legend.margin=margin(0,0,0,0),
    legend.box.margin=margin(0, 0, 20, 0),
    axis.text = element_text(size = rel(1.2)),
    axis.title = element_text(size = rel(1.35)),
    plot.caption = element_text(size = rel(1.15)))

ggsave("figures/appendix_figure_e3.pdf", 
       height = 8, width = 10, dpi = 1600)

### E.4 Reviewer-SuR Jackknife Bootstrapping -------------------------------------------------------
#### Let's do the jackknife
#### Each simulation requires more than 6 hours to complete.

sample_df <- as.data.frame(sample)

sample |> 
  dplyr::select(
    aa_lexi, da_lexi, ad_lexi,
    politysenderminustarget, 
    share_combined,
    alliance, idealdistance, reviewer_revieweed,
    to_mem, from_mem, regionalism,
    sender_hrs, target_hrs,
    migrantscat, socioecon, vulnerable,
    womenchild, physint, justicecat, political, to_cow, from_cow, year) |> 
  drop_na() |> pull(to_cow) |> unique() |> sort() -> uniqueToCOW

sample |> dplyr::select(aa_lexi, da_lexi, ad_lexi,
                        politysenderminustarget, 
                        share_combined,
                        alliance, idealdistance, reviewer_revieweed,
                        to_mem, from_mem, regionalism,
                        sender_hrs, target_hrs,
                        migrantscat, socioecon, vulnerable,
                        womenchild, physint, justicecat, political, to_cow, from_cow, year) |> 
  drop_na() |> pull(from_cow) |> unique() |> sort() -> uniqueFromCOW
pb <- progress_estimated(length(uniqueToCOW))
jk_ologit_model1_ToCOW <- data.frame()
for (i in 1:length(uniqueToCOW)) {
  pb$tick()$print()
  tryCatch({
    model1 <- MASS::polr(
      as.factor(severity) ~ 
        aa_lexi + da_lexi + ad_lexi +
        politysenderminustarget + 
        share_combined +
        alliance + idealdistance + reviewer_revieweed +
        to_mem + from_mem + regionalism +
        sender_hrs + target_hrs +
        migrantscat + socioecon + vulnerable +
        womenchild + physint + justicecat + political +
        as.factor(to_cow) + as.factor(from_cow) +
        as.factor(year),
      data = sample_df |> dplyr::filter(!to_cow %in% uniqueToCOW[i]), Hess = T,
      method = "logistic")
    
    coef <- MASS::mvrnorm(n = 4000, mu = c(coef(model1), model1$zeta), Sigma = vcov(model1))
    temp_df <- data.frame(
      term = c("aa_lexi", "da_lexi", "ad_lexi"),
      median = c(mean(coef[, "aa_lexi"], na.rm = T), 
                 mean(coef[, "da_lexi"], na.rm = T),
                 mean(coef[, "ad_lexi"], na.rm = T)),
      ll = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.025), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.025),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.025)),
      ul = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.975), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.975),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
      mutate(Exclude = paste0("SUR ", uniqueToCOW[i], " excluded"),
             id = paste0(uniqueToCOW[i]))
    jk_ologit_model1_ToCOW <- jk_ologit_model1_ToCOW |> bind_rows(temp_df)
  }, error = function(e) {
    message(paste("Error in iteration", i, ":", e$message))
    # Optionally, you can log or handle the error here
  })
}

jk_ologit_model1_FromCOW <- data.frame()
pb <- progress_estimated(length(uniqueFromCOW))
for (i in 1:length(uniqueFromCOW)) {
  pb$tick()$print()
  tryCatch({
    model1 <- MASS::polr(
      as.factor(severity) ~ 
        aa_lexi + da_lexi + ad_lexi +
        politysenderminustarget + 
        share_combined +
        alliance + idealdistance + reviewer_revieweed +
        to_mem + from_mem + regionalism +
        migrantscat + socioecon + vulnerable +
        womenchild + physint + justicecat + political +
        as.factor(to_cow) + as.factor(from_cow) +
        as.factor(year),
      data = sample_df |> dplyr::filter(!to_cow %in% uniqueFromCOW[i]), Hess = T,
      method = "logistic")
    
    coef <- MASS::mvrnorm(n = 4000, mu = c(coef(model1), model1$zeta), Sigma = vcov(model1))
    temp_df <- data.frame(
      term = c("aa_lexi", "da_lexi", "ad_lexi"),
      median = c(mean(coef[, "aa_lexi"], na.rm = T), 
                 mean(coef[, "da_lexi"], na.rm = T),
                 mean(coef[, "ad_lexi"], na.rm = T)),
      ll = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.025), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.025),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.025)),
      ul = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.975), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.975),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
      mutate(Exclude = paste0("Reviewer ", uniqueFromCOW[i], " excluded"),
             id = paste0(uniqueFromCOW[i]))
    jk_ologit_model1_FromCOW <- jk_ologit_model1_FromCOW |> bind_rows(temp_df)
  }, error = function(e) {
    message(paste("Error in iteration", i, ":", e$message))
    # Optionally, you can log or handle the error here
  })
}

### Trade lag -1 ----------
gc()
sample |> 
  dplyr::select(aa_lexi, da_lexi, ad_lexi,
                politysenderminustarget, 
                lagtrade1,alliance, idealdistance, reviewer_revieweed,
                to_mem, from_mem, regionalism,
                sender_hrs, target_hrs,
                migrantscat, socioecon, vulnerable,
                womenchild, physint, justicecat, political, to_cow, from_cow, year) |> 
  drop_na() |> pull(to_cow) |> unique() |> sort() -> uniqueToCOWt1

sample |> dplyr::select(aa_lexi, da_lexi, ad_lexi,
                        politysenderminustarget, 
                        lagtrade1,
                        alliance, idealdistance, reviewer_revieweed,
                        to_mem, from_mem, regionalism,
                        sender_hrs, target_hrs,
                        migrantscat, socioecon, vulnerable,
                        womenchild, physint, justicecat, political, to_cow, from_cow, year) |> 
  drop_na() |> pull(from_cow) |> unique() |> sort() -> uniqueFromCOWt1

jk_ologit_model3_ToCOW <- data.frame()
pb <- progress_estimated(length(uniqueToCOWt1))

for (i in 1:length(uniqueToCOWt1)) {
  pb$tick()$print()
  tryCatch({
    model1 <- MASS::polr(
      as.factor(severity) ~ 
        aa_lexi + da_lexi + ad_lexi +
        politysenderminustarget + 
        lagtrade1 +
        alliance + idealdistance + reviewer_revieweed +
        to_mem + from_mem + regionalism +
        migrantscat + socioecon + vulnerable +
        womenchild + physint + justicecat + political +
        as.factor(to_cow) + as.factor(from_cow) +
        as.factor(year),
      data = sample_df |> dplyr::filter(!to_cow %in% uniqueToCOWt1[i]), Hess = T,
      method = "logistic")
    
    coef <- MASS::mvrnorm(n = 4000, mu = c(coef(model1), model1$zeta), Sigma = vcov(model1))
    temp_df <- data.frame(
      term = c("aa_lexi", "da_lexi", "ad_lexi"),
      median = c(mean(coef[, "aa_lexi"], na.rm = T), 
                 mean(coef[, "da_lexi"], na.rm = T),
                 mean(coef[, "ad_lexi"], na.rm = T)),
      ll = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.025), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.025),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.025)),
      ul = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.975), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.975),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
      mutate(Exclude = paste0("SUR ", uniqueToCOWt1[i], " excluded"),
             id = paste0(uniqueToCOWt1[i]))
    jk_ologit_model3_ToCOW <- jk_ologit_model3_ToCOW |> bind_rows(temp_df)
    
  }, error = function(e) {
    message(paste("Error in iteration", i, ":", e$message))
    # Optionally, you can log or handle the error here
  })
}

jk_ologit_model3_FromCOW <- data.frame()
pb <- progress_estimated(length(uniqueFromCOWt1))

for (i in 1:length(uniqueFromCOWt1)) {
  pb$tick()$print()
  tryCatch({
    model1 <- MASS::polr(
      as.factor(severity) ~ 
        aa_lexi + da_lexi + ad_lexi +
        politysenderminustarget + 
        lagtrade3 +
        alliance + idealdistance + reviewer_revieweed +
        to_mem + from_mem + regionalism +
        migrantscat + socioecon + vulnerable +
        womenchild + physint + justicecat + political +
        as.factor(to_cow) + as.factor(from_cow) +
        as.factor(year),
      data = sample_df |> dplyr::filter(!to_cow %in% uniqueFromCOWt1[i]), Hess = T,
      method = "logistic")
    
    coef <- MASS::mvrnorm(n = 4000, mu = c(coef(model1), model1$zeta), Sigma = vcov(model1))
    temp_df <- data.frame(
      term = c("aa_lexi", "da_lexi", "ad_lexi"),
      median = c(mean(coef[, "aa_lexi"], na.rm = T), 
                 mean(coef[, "da_lexi"], na.rm = T),
                 mean(coef[, "ad_lexi"], na.rm = T)),
      ll = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.025), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.025),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.025)),
      ul = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.975), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.975),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
      mutate(Exclude = paste0("Reviewer ", uniqueFromCOWt1[i], " excluded"),
             id = paste0(uniqueFromCOWt1[i]))
    jk_ologit_model3_FromCOW <- jk_ologit_model3_FromCOW |> bind_rows(temp_df)
    
  }, error = function(e) {
    message(paste("Error in iteration", i, ":", e$message))
    # Optionally, you can log or handle the error here
  })
}
### Trade lag -3 ----------

sample |> dplyr::select(aa_lexi, da_lexi, ad_lexi,
                        politysenderminustarget, 
                        lagtrade3,
                        alliance, idealdistance, reviewer_revieweed,
                        to_mem, from_mem, regionalism,
                        sender_hrs, target_hrs,
                        migrantscat, socioecon, vulnerable,
                        womenchild, physint, justicecat, political, to_cow, from_cow, year) |> 
  drop_na() |> pull(to_cow) |> unique() |> sort() -> uniqueToCOWt3

sample |> dplyr::select(aa_lexi, da_lexi, ad_lexi,
                        politysenderminustarget, 
                        lagtrade3,
                        alliance, idealdistance, reviewer_revieweed,
                        to_mem, from_mem, regionalism,
                        sender_hrs, target_hrs,
                        migrantscat, socioecon, vulnerable,
                        womenchild, physint, justicecat, political, to_cow, from_cow, year) |> 
  drop_na() |> pull(from_cow) |> unique() |> sort() -> uniqueFromCOWt3

jk_ologit_model4_ToCOW <- data.frame()
pb <- progress_estimated(length(uniqueToCOWt3))

for (i in 1:length(uniqueToCOWt3)) {
  pb$tick()$print()
  tryCatch({
    model1 <- MASS::polr(
      as.factor(severity) ~ 
        aa_lexi + da_lexi + ad_lexi +
        politysenderminustarget + 
        lagtrade3 +
        alliance + idealdistance + reviewer_revieweed +
        to_mem + from_mem + regionalism +
        migrantscat + socioecon + vulnerable +
        womenchild + physint + justicecat + political +
        as.factor(to_cow) + as.factor(from_cow) +
        as.factor(year),
      data = sample_df |> dplyr::filter(!to_cow %in% uniqueToCOWt3[i]), Hess = T,
      method = "logistic")
    
    coef <- MASS::mvrnorm(n = 4000, mu = c(coef(model1), model1$zeta), Sigma = vcov(model1))
    temp_df <- data.frame(
      term = c("aa_lexi", "da_lexi", "ad_lexi"),
      median = c(mean(coef[, "aa_lexi"], na.rm = T), 
                 mean(coef[, "da_lexi"], na.rm = T),
                 mean(coef[, "ad_lexi"], na.rm = T)),
      ll = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.025), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.025),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.025)),
      ul = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.975), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.975),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
      mutate(Exclude = paste0("SUR ", uniqueToCOWt3[i], " excluded"),
             id = paste0(uniqueToCOWt3[i]))
    jk_ologit_model4_ToCOW <- jk_ologit_model4_ToCOW |> bind_rows(temp_df)
    
  }, error = function(e) {
    message(paste("Error in iteration", i, ":", e$message))
    # Optionally, you can log or handle the error here
  })
}

jk_ologit_model4_FromCOW <- data.frame()
pb <- progress_estimated(length(uniqueFromCOWt3))

for (i in 1:length(uniqueFromCOWt3)) {
  pb$tick()$print()
  tryCatch({
    model1 <- MASS::polr(
      as.factor(severity) ~ 
        aa_lexi + da_lexi + ad_lexi +
        politysenderminustarget + 
        lagtrade3 +
        alliance + idealdistance + reviewer_revieweed +
        to_mem + from_mem + regionalism +
        migrantscat + socioecon + vulnerable +
        womenchild + physint + justicecat + political +
        as.factor(to_cow) + as.factor(from_cow) +
        as.factor(year),
      data = sample_df |> dplyr::filter(!to_cow %in% uniqueFromCOWt3[i]), Hess = T,
      method = "logistic")
    
    coef <- MASS::mvrnorm(n = 4000, mu = c(coef(model1), model1$zeta), Sigma = vcov(model1))
    temp_df <- data.frame(
      term = c("aa_lexi", "da_lexi", "ad_lexi"),
      median = c(mean(coef[, "aa_lexi"], na.rm = T), 
                 mean(coef[, "da_lexi"], na.rm = T),
                 mean(coef[, "ad_lexi"], na.rm = T)),
      ll = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.025), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.025),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.025)),
      ul = c(quantile(coef[, "aa_lexi"], na.rm = T, probs = 0.975), 
             quantile(coef[, "da_lexi"], na.rm = T, probs = 0.975),
             quantile(coef[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
      mutate(Exclude = paste0("Reviewer ", uniqueFromCOWt3[i], " excluded"),
             id = paste0(uniqueFromCOWt3[i]))
    jk_ologit_model4_FromCOW <- jk_ologit_model4_FromCOW |> bind_rows(temp_df)
    
  }, error = function(e) {
    message(paste("Error in iteration", i, ":", e$message))
    # Optionally, you can log or handle the error here
  })
}

save(jk_ologit_model1_FromCOW, jk_ologit_model1_ToCOW, 
     file = "data/estimation/jk_ologit_sims_base.rdata")

save(jk_ologit_model3_FromCOW, jk_ologit_model3_ToCOW,
     file = "data/estimation/jk_ologit_sims_tradet1.rdata")

save(jk_ologit_model4_FromCOW, jk_ologit_model4_ToCOW,
     file = "data/estimation//jk_ologit_sims_tradet3.rdata")

coef_ologit1 <- MASS::mvrnorm(n = 4000, 
                              mu = c(coef(ologit_model2), ologit_model2$zeta), 
                              Sigma = vcov(ologit_model2))
ologit1 <- data.frame(
  term = c("aa_lexi", "da_lexi", "ad_lexi"),
  median = c(mean(coef_ologit1[, "aa_lexi"], na.rm = T), 
             mean(coef_ologit1[, "da_lexi"], na.rm = T),
             mean(coef_ologit1[, "ad_lexi"], na.rm = T)),
  ll = c(quantile(coef_ologit1[, "aa_lexi"], na.rm = T, probs = 0.025), 
         quantile(coef_ologit1[, "da_lexi"], na.rm = T, probs = 0.025),
         quantile(coef_ologit1[, "ad_lexi"], na.rm = T, probs = 0.025)),
  ul = c(quantile(coef_ologit1[, "aa_lexi"], na.rm = T, probs = 0.975), 
         quantile(coef_ologit1[, "da_lexi"], na.rm = T, probs = 0.975),
         quantile(coef_ologit1[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
  mutate(Exclude = paste0("No Reviewer excluded"),
         id = "main")

coef_ologit2 <- MASS::mvrnorm(n = 4000, 
                              mu = c(coef(ologit_model3), ologit_model3$zeta), 
                              Sigma = vcov(ologit_model3))
ologit2 <- data.frame(
  term = c("aa_lexi", "da_lexi", "ad_lexi"),
  median = c(mean(coef_ologit2[, "aa_lexi"], na.rm = T), 
             mean(coef_ologit2[, "da_lexi"], na.rm = T),
             mean(coef_ologit2[, "ad_lexi"], na.rm = T)),
  ll = c(quantile(coef_ologit2[, "aa_lexi"], na.rm = T, probs = 0.025), 
         quantile(coef_ologit2[, "da_lexi"], na.rm = T, probs = 0.025),
         quantile(coef_ologit2[, "ad_lexi"], na.rm = T, probs = 0.025)),
  ul = c(quantile(coef_ologit2[, "aa_lexi"], na.rm = T, probs = 0.975), 
         quantile(coef_ologit2[, "da_lexi"], na.rm = T, probs = 0.975),
         quantile(coef_ologit2[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
  mutate(Exclude = paste0("No SUR excluded"),
         id = "main1")


coef_ologit3 <- MASS::mvrnorm(n = 4000, 
                              mu = c(coef(ologit_model4), ologit_model4$zeta), 
                              Sigma = vcov(ologit_model4))
ologit3 <- data.frame(
  term = c("aa_lexi", "da_lexi", "ad_lexi"),
  median = c(mean(coef_ologit3[, "aa_lexi"], na.rm = T), 
             mean(coef_ologit3[, "da_lexi"], na.rm = T),
             mean(coef_ologit3[, "ad_lexi"], na.rm = T)),
  ll = c(quantile(coef_ologit3[, "aa_lexi"], na.rm = T, probs = 0.025), 
         quantile(coef_ologit3[, "da_lexi"], na.rm = T, probs = 0.025),
         quantile(coef_ologit3[, "ad_lexi"], na.rm = T, probs = 0.025)),
  ul = c(quantile(coef_ologit3[, "aa_lexi"], na.rm = T, probs = 0.975), 
         quantile(coef_ologit3[, "da_lexi"], na.rm = T, probs = 0.975),
         quantile(coef_ologit3[, "ad_lexi"], na.rm = T, probs = 0.975))) |> 
  mutate(Exclude = paste0("No SUR excluded"),
         id = "main2")

bind_rows(
  ologit1 |> mutate(model = "Benchmark\nModel"),
  ologit2 |> mutate(model = "Model\nwith\nTrade 1Y Lag"),
  ologit3 |> mutate(model = "Model\nwith\nTrade 3Y Lag"),
  jk_ologit_model1_FromCOW |> mutate(model = "Benchmark\nModel", jk = "Reviewer Jackknife"),
  jk_ologit_model1_ToCOW |> mutate(model = "Benchmark\nModel", jk = "SUR Jackknife"),
  jk_ologit_model3_FromCOW |> mutate(model = "Model\nwith\nTrade 1Y Lag", jk = "Reviewer Jackknife"),
  jk_ologit_model3_ToCOW |> mutate(model = "Model\nwith\nTrade 1Y Lag", jk = "SUR Jackknife"),
  jk_ologit_model4_FromCOW |> mutate(model = "Model\nwith\nTrade 3Y Lag", jk = "Reviewer Jackknife"),
  jk_ologit_model4_ToCOW |> mutate(model = "Model\nwith\nTrade 3Y Lag", jk = "SUR Jackknife")
) |> group_by(model, jk, term) |> summarize(mean = mean(median, na.rm = T)) |> print(n = Inf)

bind_rows(
  ologit1 |> mutate(model = "Benchmark\nModel"),
  ologit2 |> mutate(model = "Model\nwith\nTrade 1Y Lag"),
  ologit3 |> mutate(model = "Model\nwith\nTrade 3Y Lag"),
  jk_ologit_model1_FromCOW |> mutate(model = "Benchmark\nModel", jk = "Reviewer Jackknife"),
  jk_ologit_model1_ToCOW |> mutate(model = "Benchmark\nModel", jk = "SUR Jackknife"),
  jk_ologit_model3_FromCOW |> mutate(model = "Model\nwith\nTrade 1Y Lag", jk = "Reviewer Jackknife"),
  jk_ologit_model3_ToCOW |> mutate(model = "Model\nwith\nTrade 1Y Lag", jk = "SUR Jackknife"),
  jk_ologit_model4_FromCOW |> mutate(model = "Model\nwith\nTrade 3Y Lag", jk = "Reviewer Jackknife"),
  jk_ologit_model4_ToCOW |> mutate(model = "Model\nwith\nTrade 3Y Lag", jk = "SUR Jackknife")
) |> mutate(
  term = factor(term,
                levels = c("aa_lexi", "ad_lexi", "da_lexi"),
                labels = c("Autocratic Dyads", 
                           "Auto(REV)-Demo(SUR) Mixed Dyads",
                           "Demo(REV)-Auto(SUR) Mixed Dyads"))) -> 
  full_set


bind_rows(
  ologit1 |> mutate(model = "Benchmark\nModel"),
  ologit2 |> mutate(model = "Model\nwith\nTrade 1Y Lag"),
  ologit3 |> mutate(model = "Model\nwith\nTrade 3Y Lag")) |> 
  mutate(
    term = factor(
      term,
      levels = c("aa_lexi", "ad_lexi", "da_lexi"),
      labels = c("Autocratic Dyads",  
                 "Auto(REV)-Demo(SUR) Mixed Dyads",
                 "Demo(REV)-Auto(SUR) Mixed Dyads"))) -> 
  main_results_tab

main_results_tab |> mutate(upper = ul + 0.12) -> main_results_tab

full_set |> 
  ggplot(aes(x = model, y = median)) +
  geom_linerange(aes(ymin = ll, ymax = ul, color = id), position = position_dodge2(.5), show.legend = F) +
  geom_pointrange(aes(ymin = ll, ymax = ul, color = id), position = position_dodge2(.5), show.legend = F) +
  geom_point(aes(y = median, color = id), position = position_dodge2(.5), linetype = "dashed", show.legend = F) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey20") +
  gghighlight::gghighlight(id %in% c("main", "main1", "main2"),
                           calculate_per_facet = T,
                           use_direct_label = F,
                           unhighlighted_params = list(color = "grey80", fill = "white")) +
  geom_text(data = main_results_tab, 
            aes(x = model, y = upper, label = round(median, 2), color = id), 
            position = position_dodge2(.5),
            show.legend = F) +
  labs(y = "Estimates\n", x = "\n",
       caption = str_wrap("Note: The colored estimates for key explanatory variables are derived from models of benchmark, with one-year lagged bilateral trades, and with three-year lagged bilateral trades. The Reviewers-SuR Jackknife analysis yielded grey-colored estimates. ", 
                          145,
                          exdent = 10)) +
  scale_y_continuous(limits = c(-1.4, 1.0), 
                     breaks = c(seq(-1.4, 1.0, 0.2)),
                     labels = scales::label_number()) + 
  scale_color_manual(values = c(futurevisions::futurevisions("mars")[1],
                                futurevisions::futurevisions("mars")[3],
                                futurevisions::futurevisions("mars")[4])) +
  facet_wrap(~term, ncol = 3) +
  theme(plot.caption = element_text(hjust = 0, size = 10))

#### Figure E.4: Estimated coefficients for the effects of different regime dyads on ---------------
####             recommendation specificity

ggsave("figures/appendix_figure_e4.pdf", height = 4, width = 10, dpi = 1600)
